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Abstract 

This paper presents an algorithm for solving the multi-dimensional 
unsteady Navier-Stokes equations for compressible flows. It is based on a 
diagonally-dominant approximate factorization procedure. The 
factorization error and the timewise linearization error associated with this 
procedure are reduced by performing Newton-type inner iterations at each 
time step. The inviscid fluxes are evaluated by the fourth-order central 
differencing scheme amended with a numerical dissipation directly 
proportional to the entire dissipative part of the truncation error intrinsic 
to the third-order-biased upwind scheme. The important features of the 
proposed solution algorithm and the finite-difference scheme are 
elucidated by the numerical results of the convection of a vortex and the 
backward-facing step flows. i 
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1. Introduction 


In the past, both upwind and central differencing schemes have been 
used in the solution algorithms for the Navier-Stokes equations. These two 
approaches have their relative merits and shortcomings. One of the 
objectives of the present effort is to find a scheme that combines the 
merits of both central differencing and upwind differencing schemes. 
Another objective is to construct a multi-dimensional solution algorithm 
which is not only time-accurate but also robust with respect to this new 
differencing scheme. 

In this work, the development of such an algorithm is based upon a 
Diagonal ly-Dominant Alternating-Direction-Implicit (DDADI) approximate 
factorization procedure [1] in conjunction with a Newton-type iterative 
process [2] for implicitly advancing the solution from one time level to the 
next. A description of this solution algorithm is given in Section 2. In 
Section 3, it is shown that the proposed difference scheme can be 
characterized as a fourth-order central differencing scheme with its added 
dissipation being consistent with the numerical dissipation intrinsic to the 
third-order-biased upwind scheme, and is termed here as the FCTD 
differencing scheme (Fourth-order Central with Truncation-error-type 
Dissipation). 

The convection of an inviscid vortex in free stream has been used to 
study the performances of the proposed algorithm/scheme. Some selected 
results of this test case are presented in Section 4. To demonstrate its 
applicability for the viscous flows, the proposed algorithm/scheme has also 
been used to calculate the flows over a backward-facing step. The 
comparisons between the calculated results and the corresponding 
experimental data are shown in Section 5. The stability characteristics of 
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the proposed algorithm/scheme when applied to a three-dimensional 
domain are described in the Appendix. 

2. The Iterative Implicit Diagonally-Dominant 
Approximate Factorization 

This is the overall algorithm for solving the unsteady multi-dimensional 
compressible Navier-Stokes equations. To describe this algorithm, we 
consider the two-dimensional Navier-Stokes equations written in 
generalized coordinates (^, t|). 


^ = 0 

dr d? dn (1) 

where Q = QU ; Q ={p,pu,pv,ef j an d J = £ x Tl y -£ y ri x is the metric Jacobian. 
Here, t is the time; p is the fluid density; e is the total internal energy per 
unit volume; u and v are the velocity components in the x and y directions 
of a Cartesian coordinates system. The transformed inviscid fluxes are 

A /N ^ 

denoted by E and F\ and the transformed viscous fluxes are denoted by Ey 

A 

and Fy. The specific forms of these transformed fluxes are well known and 
will not be repeated here. 

To advance the solution of Eq. (1) from the n-th time level to the (n+1)- 
th level, the iterative implicit technique of Ref. [2] is adopted. Furthermore, 
the three-point backward time differencing is used to attain the second- 
order temporal accuracy. Let 

(Sff = /' +1) - l = 1,2,3,- m (2) 
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(3) 


(QjP = Q (n) , / = 1 

i(4 l] MsQf n) = Q {n+1) -Q' n \ 

1=1 (4) 

where f denotes an arbitrary quantity, / represents an iteration index and 
m is an intermediate iteration level between the n-th and the (n+l)-th 
time levels. The application of the iterative implicit technique to Eq. (1) 
then yields 


-CLSQ + A[gE - + -^<5F - <F„) 

^ fr] 


K m ) 


= KHS {m) 


(5) 


where a= 1.5 and 


RHS ( m ) = --1 
A t | 


+{l-atsoj n '^ 

^E-E v ) + ^F-F v f m) 

% dr] 


( 6 ) 


Eqs. (5) and (6) represent a Newton-type of timewise linearization on 

p.(n+l) 

Eq. (1) with Q as the independent variable. By iterating to 
convergence, the linearization errors associated with the residuals of the 
RHS can be driven to zero. Thus, on its convergence, the RHS approaches to 
a fully implicit nonlinear backward time differencing approximation of the 
entire unsteady Navier-Stokes equations. The numerical evaluation of the 
fluxes in Eq. (6) is discussed in Section 3. In the following, the construction 
of operators approximating the variation of fluxes appearing in the left 
hand side of Eq. (5) is described. 
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Associated with each inviscid flux vector, there exists a Jacobian matrix 
which can be further split into two parts, one is associated with 
nonnegative eigenvalues and the other with nonpositive eigenvalues. For 
example, 


A 




(7) 


are the split Jacobian matrices associated with the flux vector E. For 
convenience, the symbol ' A ' will be dropped from the flow quantities in the 
subsequent equations and discussion. 

By applying the first-order upwind split-flux technique in the context of 
a Newton-type linearization procedure, the operator approximating the 
variation of the convective flux is of the following form: 

^k 5E Xj = + [ A tj ■ A i, A^kj + A i+l,A 5 Q)i+l,j\ 

oq ( 8 ) 


In Ref. [3], the matrices A + and A' are evaluated according to 
Ai+ij = Ai + ij j A i+lj A i+lj 

and 

Aij = Aij , A i,j - A i,j (9) 

where A H and A ’ are defined by Eq. (7). In the present work, A + and A are 
evaluated at the Roe-averaged state [4]. For example, 

Aitij = A + {i - 1, i ; £) (10a) 
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which means that the flow variables are evaluated at the Roe-averaged 
state between the nodal points (i-1, j) and (i, j) while the metrics are 
evaluated at the nodal point (i, j), rather than at the intermediate point 
(i-1/2 J). The rest of the matrices are evaluated according to 


Ay = A + (t - 1, i ; §•) = A-ij 


(10b) 


Aij - A (i, i + 1, ; £•) = Ai + ij 


(10c) 


The same technique is applied to evaluate ^/dri j n w hich the 

associated split matrices are denoted by B + and B . Thus, the inviscid terms 
in Eq. (5) yield 


^SE)i,j + ii&kj = [ D + L + U](SQ)ij 
dt, drj 


(ID 


where the operators [D], [L], and [U] are defined as 


[D]( ■ 

Kj - { A tj ■ A i, j + B tj ■ B i,j)( ■ H,j 

(12a) 

[L](- 

Kj - {- A i-lj)( ■ H-l J + ■ )i,j-l 

(12b) 

[U]( ■ 

Kj - { A i+l,j)( ■ H+l J + { B i,j+l)( ■ )i,j+l 

(12c) 


The construction of the operator approximating the variation of the 
viscous fluxes follows the approach of a pointwise Jacobi-iteration 
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procedure, in which only those parts directly adding terms to the diagonal 
operator, i.e. Eq. (12a), are retained. The detailed forms of these terms are 
rather involved and are not shown here. Let V represent the Jacobian of 
these terms, then, the viscous terms in Eq. (5) yield 

-?4-8Ev) + -k-SFv) = [ S\SQ)iJ 

dt, di) ( 13 ) 

where the operator [S] is defined as 

[s]( ■ )i,j = (Vij)( ■ )i,j (14) 

It is further noted here that, if there are external source terms due to, e.g., 
gravitational acceleration or heat addition, these terms will be, along with 
Vij, included in the operator [S]. 

Let I be the identity operator, the evolution process of Eq.(5) is then 
approximated by 

[oi + MS + D+L + U)](SQ)iJ = Ar(RHS)ij ( 15) 

In practice, Eq. (15) is solved by the approximate factorization procedure. 
In addition, the diagonally-dominant treatment suggested in Ref. [1] is 
adopted. However, the splitting technique used in the present work is 
different from that suggested in Ref. [1]. A brief description of the original 
DDADI algorithm and a demonstration of its numerical instability when the 
central differencing scheme is used for the inviscid flux terms are given in 
the Appendix. In the present work, the original splitting technique is 
modified to enforce the operational symmetry of the factorization. As 
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demonstrated in the Appendix, this modification renders the overall 
solution algorithm numerically stable with respect to the central 
differencing scheme. The presently proposed splitting requires that the 
temporal iterations are carried out in pairs of advancing steps. 

Let (#?) denote the solution of the first advancing step and (#p) I+i the 
solution of the second advancing step of a pair, then (dOf at point (i,j) is 
obtained by solving 


[cd + AdS + D+L)](sQ*jj = Ax(RHS)\j 


[cd + Ad'S + D)](SQf it j = Ax(RHS)\j - [at(L + U)](sQ*f r j 

lxn ) l+1 

and V u »/ < «/ is obtained by solving 

[cd + Ai(S + D + U)](sQ* *ij = At (RHS) l fj 


[cd + AdS + D)}{SQjtj = AT(RHS) l tf - [At(U + L)](dQ* 


(16a) 

(16b) 

(17a) 

(17b) 


These equations are solved by a pointwise Jacobi-iteration procedure. 

The results obtained from several numerical experiments focused on the 
extra-long time behavior of randomly perturbed uniform background flow 
suggest that the proposed algorithm would be numerically stable, both in 
two-and three-dimensional cases, with respect to a variety of differencing 
schemes described in the following section. It is also noted here that the 
extension of this LU type of factorization to three-dimensional case is quite 
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straight forward, no additional sweep is needed to account for the added 
third dimension. 

3. The FCTD Differencing Scheme For The RHS Derivatives 
In the past, both central and upwind differencing have been employed 
to evaluate the RHS spatial derivatives. The third-order-biased upwind 
scheme (see e.g. Ref.[5]) is the lowest order scheme satisfying the 
requirement that the leading truncation error is dissipative but not 
directly contaminating the physically diffusive terms, and it has a stencil 
of five grid points in each spatial direction. The fourth-order central 
differencing also has a five-point stencil but the scheme resolution is 
higher. However, it is non-dissipative by nature, hence allows the 
excitation of spurious short-wave length modes. In order to suppress 
these spurious modes, numerical dissipation models (see e.g. [6]) containing 
fourth-difference dissipation terms are often artificially added to the finite 
difference equations. These artificial dissipation terms are not consistent 
with the Navier-Stokes equations. In the following, an analytical way of 
injecting numerical dissipation to the fourth-order central differencing 
scheme is presented. For convenience, the index j associated with the T|- 

direction is dropped in the subsequent discussion. 

Let [TU](E)i denote the third-order-biased upwind differencing of the 

first derivative of E at the spatial point (i,j), then 

[ TTJ](E)i = i-MEi+2 • Ei-i ) + L U-AEi_2 + 2 AE f.j - AEi ) 

2 h bn 


where E is the total flux, E + and E' the split fluxes, h = A£ the grid spacing 
in computational space and A is the forward differencing operator defined 

as A (‘ )* = (' )* +1 " (‘ )». Under the current effort, the flux-difference 
splitting suggested by Roe [4] is used in Eq.(18). 

A Taylor's series analysis shows that the truncation error associated 
with the third-order-biased upwind differencing has two parts: a 


dissipative part (DISS) and a dispersive part (DISP). 

[ 7T/](£) t = ® + [DISS]{E)i + [DISP]{E)i 

where 

[DISS](E)i = I ift- 1 - 2) ±h 11 ^-{E + - E-\ 
1=4, Al=2 3 I’- 

and 

[DISP](E)i = I l(4 - 2 l ' 1 )j-h 1 - 1 

1=5, Al=2 3 l! \ag* Ji 


(19) 


( 20 ) 


( 21 ) 


The dissipative part, characterized by the fourth - and higher even - 
derivatives is given by Eq. (20). The dispersive part, characterized by the 
fifth - and higher odd - derivatives is given by Eq. (21). It is further noted 
here that the dissipative part depends on the derivatives of the absolute 
flux defined as I = & E , while the dispersive part depends on the 
derivatives of the total flux. 

Let [FC](E)i denote the fourth-order central differencing of the first 
derivative of E at the spatial point (i,j), then 


[«?](** - f H £ M - Ei-1 ) - ±UEi + 2 - Ei. 2 ) 

on 12 a 


( 22 ) 


It can be shown that 
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[FC](E)i = ffl + [ DISP ]* (E)i 
\oqii 


( 23 ) 


where 


[DISP)* (E)i 


£ 1(4 - 2 l ' 1 )±h 1 ' 1 

l=5,Al=2 3 y 'U 



(24) 


i.e., the truncation error associated with the fourth-order central 
differencing is entirely dispersive. By comparing Eq. (24) with Eq. (21), it 
can be seen that [DISP]*(E)i is exactly the same as [DISP](E)i. This fact, 
then, allows the establishment of a working formula for evaluating 
[DISS](E)i. 

Subtracting Eq. (23) from Eq. (19) and using the fact that 

[ DISP]* (E)t = [ DISP ] (E)i (25) 


the dissipative part of the truncation error is then given by 

[ DISS] ( E\ = [ TU] (E)i - [ FC] (E)i (26) 

Now that a closed form of [DISS](E)i becomes available, it can be used to 
inject numerical dissipation into higher-order central differencing schemes. 
Here, the scheme under concern is of the fourth-order, and the amended 
scheme is termed as the fourth-order central differencing with truncation- 
error- type dissipation (FCTD) scheme. More specifically. 


[FCTD](E)i = [FC](£); + p[DISS](E)i 


(27) 


where p is an adjustable constant. Through Eq. (26), the working formula 
for FCTD scheme is then 

[ FCTD]{E\ = (l-P)[FC}{E)i + /?[ 7T/](^ (28) 

where [TU](E)i and [FC](E)i are given by Eq. (18) and Eq. (22) respectively. 
Since 


[ FCTD]{E)i = ® + p{DISS](E)i + [DISP](E)i 
v "s li 


( 29 ) 


the effective truncation error of the proposed scheme is p[DISS](E)i + 
[DISP](E)i. It can be seen from Eq. (28) that, when p=0, FCTD becomes the 
fourth-order central differencing scheme without numerical dissipation, 
and, when p = l, it becomes the third-order-biased upwind scheme. 

In a nutshell, the proposed FCTD scheme is an amended fourth-order 
central differencing scheme with its injected numerical dissipation having 
the same form as the entire dissipative part of the truncation error 
intrinsic to the third-order-biased upwind scheme. Such a dissipation is 
essentially an infinite series with its elements being the fourth - and 
higher even - derivatives of the absolute flux. The relative amount of 
added dissipation can be controlled through an adjustable parameter p. It 
is further noted here that the FCTD scheme yields finite-difference 
approximations which are consistent with the physical flux derivatives. 

In regard to the evaluation of the viscous terms in Eq. (6), these terms 
are first written in the strong conservation form. Then, half-spacing 
second-order central differencing scheme is used to evaluate the 
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derivatives, and simple average is employed to provide mid-point values 
of transport coefficients. The details of these terms are not shown here. 

4. Results - The Convection of a Free Vortex 

In references [7,8], the convection of a Lamb-type vortex in a free 
stream was used as a test case for assessing numerical schemes' ability to 
preserve and convect concentrated vortices. This test case also is used 
here to infer the dissipative, dispersive and resolution properties of the 
proposed scheme. 

In test calculations, the reference length is the vortex core radius; the 
reference flow conditions are the free stream conditions. The number of 
grid points used is 241 x 61. Figure 1 depicts the computational domain 
and the grid distribution. In the region containing the vortex path (i.e., -5 
< x < 50, -5 < y < 5), the grid spacing is uniform with Ax = Ay=0.25. At 

t=o, the vortex center is located at the origin (x=o, y=o). The initial 
condition and the subsequent unsteady boundary conditions are specified 
through the analytic solution [8] which corresponds to the vortex being 
convected with a free stream having a Mach number AL»= 0.5. The vortex 
flow rotates in the clockwise direction. The center pressure (Pc) is a 
minimum having a value of 0.84. The farfield free stream pressure is 
P^l.O. The vorticity magnitute of this vortex is shown in Fig. 2. The 
subsequent convection of this vortex is calculated with a constant time 
step At=0.05, which is used in all the test calculations. These calculations 
are terminated when the vortex center has traveled a distance of 45 core 
radii, and this requires a total of 900 time steps. It is further noted here 
that the advancement of one time step consists of two pairs of Newton- 
type inner iterations. Based on our experience, this value of the number of 
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inner iterations represents an optimal use of computer resources for the 
present problem. 

The results obtained by applying the third-order-biased upwind scheme 
are shown in Figs. 3 and 4. The pattern of vorticity magnitude contours is 
well preserved, yet there is significant dissipation of vortex strength. It is 
noted here that the formal accuracy of this scheme is of the third order, 
and the leading truncation-error term is a dissipative fourth-derivative. 

The fourth-order central differencing scheme has the same computational 
stencil as that of the third-order-biased upwind scheme, but is formally of 
higher order accuracy. The fourth-order central differencing scheme is 
intrinsically non-dissipative, and the leading truncation -error term is a 
dispersive fifth-derivative. Therefore, it can be expected that, without the 
injection of artificial dissipation, the vortex strength would be better 
preserved, but spurious short-wave-length components would emerge. 
These features are demonstrated by results shown in Fig.5 and Fig. 6. To 

suppress the excitation of spurious components, the fourth -difference 
artificial viscosity model [6] with p *=0.015 is used (hereafter, termed as 
the FCAV scheme), and the results are presented in Figs. 7 and 8. The 
excitation of phase errors is successfully suppressed during the period of 
this simulation, and the preservation of the vortex strength is just as good 
as that shown in Fig. 6. For FCTD scheme, the coefficient of the leading 
dissipative fourth-derivative term will have a numerical value of 0.015, if 
P is chosen to be 0.18. Fig. 9 indicates that the vorticity pattern obtained 
by using the FCTD scheme is well preserved. By comparing Fig. 7(b) and 
Fig. 9(b), it can be seen that, in Fig. 7(b), there exist incipient short-wave- 
length components, while the contours in Fig. 9(b) do not exhibit these 
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spurious components. The results shown in Fig. 10 are slightly more 

dissipated than those in Fig. 8. 

The proposed FCTD scheme with (3=0.18 is more accurate than the third- 
order-biased upwind scheme, but is slightly more dissipative than the 
FCAV scheme with p *=0.01 5. However, the complete absence of spurious 
short-wave-length components from the contours of vordcity magnitude 
makes the FCTD differencing a more reliable scheme than the FCAV 
scheme. The above results are for the case of 0.5. Flows with M^= 1.0 
and M^l.6 have also been calculated to compare the performances of the 
schemes in transonic and supersonic regimes. They lead to the same 
conclusions. 


5. Results - The Backward-Facing Step Flows 
The application of the proposed algorithm/scheme to the solution of the 
full Navier-Stokes equations is demonstrated by calculating flows with 
regions of separation behind a backward-facing step mounted in a two- 
dimensional channel. The calculated results are then compared with the 
corresponding experimental data reported in Ref. [9]. This experimental 
study concluded that, in the laminar range (Re<1200), the flow will 
maintain its two-dimensionality only when Reynolds number Re<400. The 
definition of the Reynolds number is given by Re=u(2h)/v, where u is the 
average inlet velocity, which corresponds to two-thirds of the maximum 
inlet velocity, h is the height of the inlet channel, and v is the kinematic 
viscosity. In addition, the channel height downstream of the step is 
H=1.9423h, and the step height S=0.9423h. Fig. 11 depicts the size of the 
computational domain and the grid distribution. The number of grid 
points is 81x31. 


Under the current effort, all the calculations were carried out with 
Mi=0. 18 and Pr=0.72, where is the Mach number based on u, and Pr is 
the laminar Prandtl number. The calculations started with some guessed 
initial conditions and proceeded until a steady state was reached. The 
criterion for reaching an asymptotic steady state is that the maximum L 2 - 
norm residual is smaller than 5x1 O' 5 . 

Both the fourth-order central differencing with constant-coefficient 
fourth-difference artificial dissipation (J5*=0.5), i.e. the FCAV scheme, and 
the proposed FCTD scheme ((5=0.5) have been employed to calculate the 
flows with Reynolds numbers in the laminar range: Re=100, 389, and 

1000. The calculated streamwise velocity profiles at stations downstream 
of the step and their corresponding experimental data are shown in Fig. 
12(a) (Re=100), Fig. 12(b) (Re=389), and Fig. 12(c) (Re=1000). For Re=100 
and 389, the flows are experimentally two-dimensional, and the present 
two-dimensional numerical results agree very well with the measured 
data in terms of the reattachment lengths (shown in table I), and the 
streamwise velocity profiles (shown in Fig. 12(a) and 12(b)). Furthermore, 
the results obtained from the two numerical schemes are practically 
indistinguishable. At Re=1000, the experiment indicated that the flow 
loses its two-dimensionality. This is reflected in the apparent 
discrepancies between the experimental data and the two-dimensional 
numerical results as shown in Fig. 12(c). It is noted here that this kind of 
discrepancy also exists between the experimental data and two- 
dimensional numerical results obtained by using other schemes (see e.g. 
Ref. [9]). It can also be seen that, in the recirculating region, there are 
some differences between the results obtained from the two numerical 
schemes. Nevertheless, these differences are considered here as 



insignificant. 

These two schemes are further compared in terms of the contour plots 
of the static pressure (Fig. 13 and Fig. 14). For all the practical purposes, 
the corresponding contour plots can be considered as the same. However, 
they do reveal some subtle differences in the dissipative features of the 
two schemes. The FCAV scheme is more effective in smoothing out sharp 
pressure wiggles associated with the rapid variation of grid spacing in the 
middle of the larger channel. Although not shown here, it is less effective 
in damping out the small-amplitude high-wave-number modes of the 
velocity. 


6. Concluding Remarks 

The FCTD differencing scheme and an overall solution algorithm have 
been developed to solve the multi-dimensional unsteady Navier-Stokes 
equations for compressible flows. The basis of the solution algorithm is a 
diagonally-dominant approximate factorization procedure. The 
factorization error and the timewise linearization error associated with this 
baseline procedure are reduced by performing Newton-type inner 
iterations at each time step. The robustness of the overall algorithm is 
enhanced by carrying out the temporal iterations in pairs to enforce the 
operational symmetry of the factorization procedure. The temporal 
accuracy is increased to second-order by using the three-point backward 
time differencing. The viscous fluxes are evaluated by using the half- 
spacing second-order central differencing scheme. The inviscid fluxes are 
evaluated by the proposed FCTD scheme, which is an amended fourth- 
order central differencing scheme with its injected numerical dissipation 
having the same form as the entire dissipative part of the truncation-error 
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intrinsic to the third-order-biased upwind scheme. The convection of a 
free vortex and the backward-facing step flows are chosen as test cases to 
demonstrate the current capability of this solution algorithm and the FCTD 
scheme. These numerical results compare well with corresponding 
analytical solution or experimental data. 


Appendix: The DDADI Algorithm and an Investigation of Its 
Numerical Stability 

In terms of the notations used in Eq. (16) and Eq. (17), the DDADI 
algorithm suggested in Ref. [ 1 ] can be written as 


[oi + At(S + D+L)]{SQ)* = Az(RHS) 1 

[cd + 4t('S + D+[/;](5Q/ = At(RHS) 1 - OtL](aQ)* 

and 

[or + At(S + D + U)]{8Q)** = At(RHS) l+1 

[cd + Ar(S + D + L)](SQf +1 = Az(RHS) l+1 - [ArU]{8Q)~ 


(Al.a) 

(Al.b) 


(A2.a) 

(A2.b) 


where Eq. (Al) is the counterpart of Eq. (16) and Eq. (A2) is the 
counterpart of Eq. (17). For convenience, the spatial indices i,j and k have 
been dropped from these equations. By comparing the respective 

counterparts, it can be seen that the main difference lies in the treatment 
of the off-diagonal contributions at the corrector's stage. 

A measure of the factorization errors associated with a pair of 
temporal iterations described by Eqs. (Al) and (A2) is 


£ = {Azf{[L][(d + Ax(S + Z>;]' i [l/])(5Q) / 



+ (4t)*{[l/|[<tf + Ar(S + D)]- 1 [L])(dQf +1 


(A3) 


For the presently proposed algorithm described by Eqs. (16) and(17), a 
measure of the factorization errors is 


£ * = -(Axf[U] 


{ dx 


,* \ 


+ (Arf{[L][aI + Ax(S + 


-(Azf[L] 




dr 


+ (Axf {[C/][ cd + Ar(S + Z))]'^[L]}(5Q)* 


(A4) 


Upon the convergence of the temporal iterations, both e and e* are 
asymptotically removed. 

Numerical experiments were conducted to examine the stability 
characteristics of the original DDADI algorithm and the presently proposed 
algorithm. In the following, the stability characteristics of these two 
algorithms when applied to a three-dimensional domain are presented. 
Fig. A-l shows the stretched-grid layout of this rectangular domain in a 
physical space occupied by an inviscid uniform flow with Mach number 
being 0.35. When the Euler equations are numerically solved to simulate 
the evolution of this flow, initially broad-band and infinitesimal 
disturbances originating from the machine round-off error as well as the 
truncation error associated with the grid stretching will be introduced into 
the uniform background flow. The subsequent growth of these 
disturbances is closely related to the numerical stability property of the 
solution algorithm and differencing scheme. The observed long-time 



growth behavior of the disturbances serves to indicate the stability 
characteristics of the overall solution procedure. 

The domain of a 2x1x1 box shown in Fig. A-l is covered by 21x11x7 
grid points, and the background flow is in the x (i.e. Indirection. The 
boundary conditions for this investigation are 

1=1 : u = l, v = w = 0, T = l, !OL = o 

dn 2 


J = 1, K = 1 : symmetry conditions 


I = 21, J = 11, K = 7 : 


du 

dn 


du_ 

dn 


dw 

dn 


cF 

dn 


= 0, p = 1 


(A5) 


where n denotes the normal direction. 

First-order time-accurate calculations without employing the inner 
Newton-type iteration process have been conducted. The RHS inviscid flux 
terms are evaluated with the fourth-order central differencing without 
any added numerical dissipation. The initial condition is an uniform flow. 
The time-step used is a constant and has a CFL number of 250. The mid- 
point, i.e., the point at (1=11, J=6, K=4), is selected as the representative 
point for illustrating the stability characteristics. 

In the case of the original DDADI algorithm [1], the temporal growth of 
the initially infinitesimal velocity disturbances is depicted in Fig. A-2. It 
can be seen that, within 2100 time-steps, the disturbances tend to grow 
out of bounds, and this algorithm is considered as unstable with respect to 
the fourth-order central differencing scheme under the present test 
condition. 
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The solutions at the 2000-th time-step are then used as the initial 
condition for investigating the stability properties of the presently 
proposed algorithm. Part of the results are shown in Fig. A-3, which 
covers a period of the first 4000 time-steps. An additional 16000 time 
steps have been advanced. Although the entire history of convergence is 
not shown here, it is clear that converged solutions with the magnitudes of 
v-component and w-component of the order of 10 11 are obtained (machine 
accuracy is of the order of 10 ^). The above result suggests that the 
presently proposed algorithm is numerically stable in three-dimensional 
domain. 
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Fig. 1 Computational domain and grid distribution for vortex 
convection tests (only every other grid lines are shown). 



Fig. 2 Contours of the vorticity magnitude (t=0). Increment level=0.0275. 




3 Contours of the vorticity magnitude (3rd - order - biased 
upwind differencing). Increment level=0.0275. 




Fig. 4 Static pressure profile along the centerline at different 
time stations (3rd - order - biased upwind differencing). 



(a) t = 25.0 (b) t _ 45 0 


Fig. 5 Contours of the vonicity magnitude (4th - order central 

differencing without artificial dissipation). Increment levcl=0.0275. 




Fig. 6 Static pressure profile along the centerline at different 
time stations (4th - order central differencing without 
artificial dissipation). 
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(a) t = 25.0 (b) t = 45.0 

Fig. 9 Contours of the vorticity magnitude (FCTD scheme with p = 
0.18). Increment level=0.0275 
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X 


Fig. 10 Sialic pressure profile along the centerline at different 
time stations (FCTD scheme with p = 0.18). 


x=-3 



x=0 x=37 

Fig. 11 Computational domain and grid distribution for backward- 
facing step flows (only every other grid lines are 
shown). 


Re 

\msmm 

R1 

S2 

R2 

100 

msssM 

2.98 

- 


FCAV 

2.94 

- 

_ 


gxp 

3.0 

_ 


389 

FCTD 

7.71 

- 

. 

FCAV 

7.60 

- 

. 


Exp 

7.8 



500 

FCTD 

NA 

N.A 

NA 

FCAV 

8.85 

7.70 

11.70 


r ^ 

10.0 

8.2 

13.5 

600 

FCTD 

N.A 

NA 

N.A 

FCAV 

9.73 

8.04 

14.10 


ExP 

11.2 

8.6 

14.8 


FCTD 

11.72 

9.23 

22.43 

1000 

FCAV 

12.06 

9.46 

21.81 | 


Exp 

16.2 

13.4 

21.8 



Table I Measured and computed detachment and reattachment length. 
The flow pattern shown is the contours of streamwise 
velocity when Re = 1000. 
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(c) Rc = 1000 

F.g. 14 Comours of the static pressure (FCAV scheme with p* = 0.5). 



Computational domain and grid distribution for test of 
numerical stability. 






-2 Temporal evolution of the velocity disturbance at the 
mid-point (original DDADI algorithm) 


Fig. A -3 Temporal evolution of the velocity disturbance at the 
mid-point (lime-step: 1 - 4000) 
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